function ecog = ecogFindPeak(ecog, band,signaltype);

% INPUT: ecog structure and the number of the band according to the order
%        in ecog.BPdata. 
%        theta      ==> 1;
%        lowalpha   ==> 2;
%        hialpha    ==> 3;
%        beta       ==> 4;
%        gamma      ==> 5;
%        higamma    ==> 6;
%
% OUTPUT: ecog structure with a new field called peaks, which is a cell 
%         array containing the time points of the peaks found in the 
%         frequency bands specified by band. 
% 
% EXAMPLE: ecog = ecogFindpeak(ecog,1,'trial')
%          for theta peaks
%          
% NOTE: use ecogBandpassFilter before finding the peaks.  
%       Afterwards ecogValidpeaks.m can be used for detection of peaks 
%       after a certain stimulus or event

if isfield(ecog,'peaks') == 1;
    ecog = rmfield(ecog,'peaks');
end


clc;



ecog.peaks{band(1)} = zeros(size(ecog.data,1),...
                            1,...
                            size(ecog.BPdata{band(1)},3));
switch signaltype;
    case 'timeseries'
        for k = 1:size(band,2)
        for chan = 1:size(ecog.data,1);  
            sig = ecog.BPdata{1}(chan,:);
            hilsig = hilbert(sig);
            angsig = angle(hilsig);
            maxvals = find(abs(angsig) < .02);
            d = diff(maxvals);
            idy = find(d > 1);
            maxvals = maxvals(idy);
            ecog.peaks{band}(chan,1:length(maxvals)) = maxvals;

        %     d = sign(diff(ecog.BPdata{band}(chan,:)));
        %     maxvals = find(d < [d(1) d(1:end-1)]); 
        %     ecog.peaks{band}(chan,1:length(maxvals)) = maxvals;
        end
        end

    case 'trial'
        for k = 1:size(band,2);
        for trial = 1:size(ecog.BPdata{band(k)},3);
        for chan = 1:size(ecog.data,1);
            sig = ecog.BPdata{band(k)}(chan,:,trial);
%             hilsig = hilbert(sig);        
%             angsig = angle(hilsig);
%             maxvals = find(abs(angsig) < .001);          
%             d = diff(maxvals);
%             idy = find(d > 1);
%             maxvals = maxvals(idy);
    
            d = sign(diff(sig));
            maxvals = find(d < [d(1) d(1:end-1)]);
%             x = find(sig(maxvals) > mean(sig) + std(sig));
%             maxvals = maxvals(x);
            ecog.peaks{band(k)}(chan,1:length(maxvals),trial) = maxvals;
        end
        end
        end  
end